#include <interpolate.hpp>

template<uint Dim,uint N>
typename Interpolate<Dim,N>::mesh_ptrtype
Interpolate<Dim,N>::createMesh( double meshSize )
{
    // instantiate a new mesh
    mesh_ptrtype mesh( new mesh_type );

    // generate a tensorized domain (hyper cube of dimension Dim)
    GmshTensorizedDomain<convex_type::nDim,
                         convex_type::nOrder,
                         convex_type::nDim,
                         Simplex> td;
    td.setCharacteristicLength( meshSize );
    td.setX( std::make_pair( 0, 2 ) );
    std::string fname = td.generate( convex_type::name().c_str() );

    // importer the mesh generated by td
    ImporterGmsh<mesh_type> import( fname );
    mesh->accept( import );

    return mesh;
} // Interpolate::createMesh

//oooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo//
//oooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo//
//oooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo//

template<uint Dim, uint N>
typename Interpolate<Dim,N>::mesh_ptrtype
Interpolate<Dim,N>::createMesh1( double meshSize )
{
    mesh_ptrtype mesh( new mesh_type );

    std::ostringstream ostr;
    ostr << "Mesh.MshFileVersion = " << 2 << ";\n"
         << "h=" << meshSize << ";\n"
         << "a=0;\n"
         << "Point(1) = {-1,-1, 0, h};\n"
         << "Point(2) = { a,-1, 0, h};\n"
         << "Point(3) = { 0.5, 1, 0, h};\n"
         << "Point(4) = {-1, 1, 0, h};\n"
         << "Line(1) = {1, 2};\n"
         << "Line(2) = {2, 3};\n"
         << "Line(3) = {3, 4};\n"
         << "Line(4) = {4, 1};\n"
         << "Line Loop(1) = {1, 2, 3 ,4 };\n"
         << "Plane Surface(2) = {1};\n"
         << "Physical Line(\"Bord\") = {1,3,4};\n"
         << "Physical Line(\"Interface\") = {2};\n"
         << "Physical Surface(\"Physique1\") = {2};\n";


    Gmsh gmsh;
    std::string fname = gmsh.generate( "domain1", ostr.str()  );

    ImporterGmsh<mesh_type> import( fname );
    
    float gmsh_version = this->vm()["gmsh"].template as<float>();
    if (gmsh_version==2.0) import.setVersion( "2.0" );
    else if (gmsh_version==2.1) import.setVersion( "2.1" );

    mesh->accept( import );

    return mesh;

}

//oooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo//
//oooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo//
//oooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo//

template<uint Dim, uint N>
typename Interpolate<Dim,N>::mesh_ptrtype
Interpolate<Dim,N>::createMesh2( double meshSize )
{
    mesh_ptrtype mesh( new mesh_type );

    std::ostringstream ostr;
    ostr << "Mesh.MshFileVersion = " << 2 << ";\n"
         << "h=" << meshSize << ";\n"
         << "a=0;\n"
         << "Point(1) = { a,-1, 0, h};\n"
         << "Point(2) = { 1,-1, 0, h};\n"
         << "Point(3) = { 1, 1, 0, h};\n"
         << "Point(4) = { a, 1, 0, h};\n"
         << "Line(1) = {1, 2};\n"
         << "Line(2) = {2, 3};\n"
         << "Line(3) = {3, 4};\n"
         << "Line(4) = {4, 1};\n"
         << "Line Loop(1) = {1, 2, 3 ,4 };\n"
         << "Plane Surface(2) = {1};\n"
         << "Physical Line(\"Bord\") = {1,2,3};\n"
         << "Physical Line(\"Interface\") = {4};\n"
         << "Physical Surface(\"Physique1\") = {2};\n";


    Gmsh gmsh;
    std::string fname = gmsh.generate( "domain2", ostr.str()  );

    ImporterGmsh<mesh_type> import( fname );
    
    float gmsh_version = this->vm()["gmsh"].template as<float>();
    if (gmsh_version==2.0) import.setVersion( "2.0" );
    else if (gmsh_version==2.1) import.setVersion( "2.1" );

    mesh->accept( import );

    return mesh;

}

//oooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo//
//oooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo//
//oooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo//
#include <life/lifemesh/kdtree.hpp>
#include <outilsgeom.hpp>
template<uint Dim,uint N>
void
Interpolate<Dim,N>::run()
{

    if ( this->vm().count( "help" ) )
        {
            std::cout << this->optionsDescription() << "\n";
            return;
        }

    //-----------------------------------------------------------------------------------//

    int MyOrder=N;

    this->changeRepository( boost::format( "%1%/%2%/P%3%/h_%4%/" )
                            % this->about().appName()
                            % convex_type::name()
                            % MyOrder
                            % this->vm()["hsize"].template as<double>()
                            );

    //-----------------------------------------------------------------------------------//

    double meshSize = this->vm()["hsize"].template as<double>();
    mesh_ptrtype mesh1 = createMesh1( meshSize );
    mesh_ptrtype mesh2 = createMesh2( 0.07 );

    //-----------------------------------------------------------------------------------//

    //The function space and some associated elements(functions) are then defined
    space_ptrtype Xh1 = space_type::New( mesh1 );
    element_type u1( Xh1, "u1" );

    space_ptrtype Xh2 = space_type::New( mesh2 );
    element_type u2( Xh2, "u2" );

    //-----------------------------------------------------------------------------------//
    
       
/*    typedef ublas::vector<double> node_type;
    typedef boost::tuple<node_type, size_type, uint16_type > index_node_type;
    typedef std::vector<index_node_type> points_type;  
*/
    KDTree M_tree;
    
/*    // objet representant un noeud (ublas::vector<double> )
    typedef typename node<double>::type node_type;

    node_type p1(1);p1(0)=1.0;
    node_type p2(1);p1(0)=2.0;
    node_type p3(1);p1(0)=5.0;
    node_type p4(1);p1(0)=3.0;
    node_type p5(1);p1(0)=2.0;

    M_tree.addPoint(p1);
    M_tree.addPoint(p2);
    M_tree.addPoint(p3);
    M_tree.addPoint(p4);
    M_tree.addPoint(p5);
*/    
    // initialisation of the mesh::inverse data structure
    typename space_type::dof_type::dof_points_const_iterator it_dofpt = Xh1->dof()->dofPointBegin();
    typename space_type::dof_type::dof_points_const_iterator en_dofpt = Xh1->dof()->dofPointEnd();
  //  size_type nbpts = 0;
    for( ; it_dofpt != en_dofpt; ++it_dofpt /*, ++nbpts*/  ) {
      M_tree.addPointWithId( boost::get<0>( *it_dofpt ) , boost::get<1>( *it_dofpt ) , boost::get<2>( *it_dofpt ) ) ;
    }

    typename mesh_type::element_iterator el_it;
    typename mesh_type::element_iterator el_en;
    boost::tie( boost::tuples::ignore, el_it, el_en ) = Life::elements( *mesh1 );
    
    BoundingBox<> bb;

    KDTree::points_type boxpts;

    for( ; el_it != el_en; ++el_it ) {
      
        bb.make( el_it->G() );  // std::cout<<el_it->G();
        for (size_type k=0; k < bb.min.size(); ++k) { bb.min[k] -= 1e-10; bb.max[k] += 1e-10; }
        
        M_tree.pointsInBox( boxpts, bb.min, bb.max );
        
        std::cout<<"\n--------------------------------------------------------------\n";
        std::cout<<boxpts;
        std::cout<<"\n--------------------------------------------------------------\n";
        
        for( size_type i = 0; i < boxpts.size(); ++i ) {
    
            size_type index = boost::get<1>( boxpts[i] );
            std::cout<<index<< " ";
        }
      
    }


    


} // Interpolate::run

//#include <mykdtree.hpp>
//#include <localisation.hpp>
template<uint Dim,uint N>
void
Interpolate<Dim,N>::run2()
{

    if ( this->vm().count( "help" ) )
        {
            std::cout << this->optionsDescription() << "\n";
            return;
        }

    //-----------------------------------------------------------------------------------//

    int MyOrder=N;

    this->changeRepository( boost::format( "%1%/%2%/P%3%/h_%4%/" )
                            % this->about().appName()
                            % convex_type::name()
                            % MyOrder
                            % this->vm()["hsize"].template as<double>()
                            );

    //-----------------------------------------------------------------------------------//

    double meshSize = this->vm()["hsize"].template as<double>();
    mesh_ptrtype mesh1 = createMesh1( meshSize );
    mesh_ptrtype mesh2 = createMesh2( 0.07 );

    //-----------------------------------------------------------------------------------//

    //The function space and some associated elements(functions) are then defined
    space_ptrtype Xh1 = space_type::New( mesh1 );
    element_type u1( Xh1, "u1" );

    space_ptrtype Xh2 = space_type::New( mesh2 );
    element_type u2( Xh2, "u2" );

    //-----------------------------------------------------------------------------------//

    node_type pp(2);pp(0)=-0.5;pp(1)=0.5;
    //node_type pp(2);pp(0)=-3;pp(1)=0.23;

    bool bb;
    size_type indd;
    node_type __x_ref;
    typename mesh_type::Localisation LLL(mesh1);
    boost::tie( bb, indd, __x_ref)=LLL.searchElement(pp);
    std::cout<<"\nrecherche1="<<indd;

    typename mesh_type::Localisation LLL2(mesh1);
    boost::tie( bb, indd, __x_ref)=LLL2.searchElement(pp);
    std::cout<<"\nrecherche2="<<indd;


    KDTree kd_tree;


    typedef typename node<double>::type node_type;

    
    typename mesh_type::element_iterator el_it;
    typename mesh_type::element_iterator el_en;
    boost::tie( boost::tuples::ignore, el_it, el_en ) = Life::elements( *mesh1 );

    typedef boost::tuple<node_type, std::list<size_type> > node_elem_type;

    //std::map<size_type, std::list<size_type> > M_geoGlob_Elts;
    std::map<size_type, node_elem_type > M_geoGlob_Elts;
    
    std::cout<<"\nSizeMap : "<<M_geoGlob_Elts.size();

    for( ; el_it != el_en; ++el_it ) {
        for (int i=0;i<el_it->nPoints();++i) {
          if (boost::get<1>( M_geoGlob_Elts[el_it->point(i).id()] ).size()==0) {
              boost::get<0>( M_geoGlob_Elts[el_it->point(i).id()] ) = el_it->point(i).node();
           
           kd_tree.addPoint(el_it->point(i).node(),el_it->point(i).id() );
           }
          boost::get<1>( M_geoGlob_Elts[el_it->point(i).id()] ).push_back(el_it->id());
        }
    }

    /*
    std::map<size_type, node_elem_type >::iterator it=M_geoGlob_Elts.begin();
    std::map<size_type, node_elem_type >::iterator en=M_geoGlob_Elts.end();

    for( ; it != en; ++it ) {  std::cout<< boost::get<0>(it->second)<<"\n";   }
    */
      

    kd_tree.search(pp);
    kd_tree.showResultSearch();
    
    typename KDTree::points_search_type ptsNN = kd_tree.pointsNearNeighbor();

    typename KDTree::points_search_const_iterator itNN = ptsNN.begin();
    typename KDTree::points_search_const_iterator itNN_end = ptsNN.end();
    
    typename std::list<size_type>::iterator itL;
    typename std::list<size_type>::iterator itL_end;
    
    for( ; itNN != itNN_end; ++itNN ) {
        std::cout<<"\n"<<boost::get<3>( *itNN )<< " : ";
        itL= boost::get<1>( M_geoGlob_Elts[boost::get<3>( *itNN )] ).begin();
        itL_end= boost::get<1>( M_geoGlob_Elts[boost::get<3>( *itNN )] ).end();
        for( ; itL != itL_end; ++itL ) { std::cout<<*itL<<" "; }
        std::cout<<"\n";
    }
    
    std::list< std::pair<size_type, uint> > ListTri;
    std::list< std::pair<size_type, uint> >::iterator itLF;
    std::list< std::pair<size_type, uint> >::iterator itLF_end;
    
    itNN = ptsNN.begin();
    itNN_end = ptsNN.end();
    
    for( ; itNN != itNN_end; ++itNN ) {
        itL= boost::get<1>( M_geoGlob_Elts[boost::get<3>( *itNN )] ).begin();
        itL_end= boost::get<1>( M_geoGlob_Elts[boost::get<3>( *itNN )] ).end();
        for( ; itL != itL_end; ++itL ) {
            itLF=ListTri.begin();
            itLF_end=ListTri.end();bool find=false;
            while (itLF != itLF_end && !find ) {
                if (itLF->first == *itL) find=true;
                else ++itLF;
            }
            if (find) {
                uint nb=itLF->second+1;
                size_type numEl=itLF->first;
              //  ListTri.erase(itLF);
		    ListTri.remove(*itLF);
                itLF=ListTri.begin();
                itLF_end=ListTri.end();bool find=false;
                while (itLF != itLF_end && !find ) {
                    if (itLF->second < nb) find=true;
                    else ++itLF;
                }
                ListTri.insert(itLF,std::make_pair(numEl,nb));
            
            }
            else ListTri.push_back(std::make_pair(*itL,1));

                
        }
        
    }
    
    itLF=ListTri.begin();
    itLF_end=ListTri.end();
    std::cout<<"\nApres trie\n-----------\n";
    for( ; itLF != itLF_end; ++itLF ) {
        std::cout<<itLF->first<<" : "<<itLF->second<<"\n";
    }


    typename mesh_type::element_type elt=mesh1->element((ListTri.begin())->first );
 
   //-----------------------------------------------------------------------------------//
    
    typedef typename mesh_type::element_type element_mesh_type;
    typedef typename mesh_type::gm_type::template Context<vm::JACOBIAN|vm::KB|vm::POINT, element_mesh_type> gmc_type;
    typedef boost::shared_ptr<gmc_type> gmc_ptrtype;
    typename mesh_type::gm_type::reference_convex_type refelem;
    typename mesh_type::gm_type::precompute_ptrtype __geopc( new typename mesh_type::gm_type::precompute_type( mesh1->gm(), refelem.points() ) );

    
    
    gmc_ptrtype __c( new gmc_type( mesh1->gm(), elt, __geopc ) );
                                   
                                   
    // node_type q(2);q(0)=-0.5;q(1)=0.0;                                   
                                   
    // for( ; el_it != el_en; ++el_it ) {
        // get geometric transformation
        __c->update( elt );
        typename mesh_type::Inverse::gic_type gic( mesh1->gm(), elt );
        
        gic.setXReal( pp);//boost::get<0>( boxpts[i] ) );
        bool isin;
        double dmin;
        boost::tie( isin, dmin ) = refelem.isIn( gic.xRef() );
        
        if (isin) std::cout<<"\ntrue0\n"; 
	else std::cout<<"\nfalse0\n";




    //-----------------------------------------------------------------------------------//
    /*    
    boost::tie( boost::tuples::ignore, el_it, el_en ) = Life::elements( *mesh1 );
    
    typedef typename mesh_type::element_type element_mesh_type;
    typedef typename mesh_type::gm_type::template Context<vm::JACOBIAN|vm::KB|vm::POINT, element_mesh_type> gmc_type;
    typedef boost::shared_ptr<gmc_type> gmc_ptrtype;
    typename mesh_type::gm_type::reference_convex_type refelem;
    typename mesh_type::gm_type::precompute_ptrtype __geopc( new typename mesh_type::gm_type::precompute_type( mesh1->gm(), refelem.points() ) );

    
    
    gmc_ptrtype __c( new gmc_type( mesh1->gm(), *el_it, __geopc ) );
                                   
                                   
    node_type q(2);q(0)=-0.5;q(1)=0.0;                                   
                                   
    for( ; el_it != el_en; ++el_it ) {
        // get geometric transformation
        __c->update( *el_it );
        typename mesh_type::Inverse::gic_type gic( mesh1->gm(), *el_it );
        
        gic.setXReal( q);//boost::get<0>( boxpts[i] ) );
        bool isin;
        double dmin;
        boost::tie( isin, dmin ) = refelem.isIn( gic.xRef() );
        
//        if (isin) std::cout<<"\ntrue\n";
//        else std::cout<<"\nfalse\n";
    }                             
    */                             
/*
    std::cout<<"\n\nnumPoints : "<< mesh1->numPoints();

    std::cout<<"\n\nnumPoints : "<< mesh1->point(0).node();

    std::cout<<"\n\nnumPoints/Elem : "<< el_it->nPoints();
    
    std::cout<<"\n\nnumPoints/Elem : "<< el_it->point(0).id();
    
 for (int i=1;i<mesh1->numPoints();++i)
     std::cout<<"\n\nnumPoints : "<< mesh1->point(i).id();
  */  

/*
    typename space_type::dof_type::dof_points_const_iterator it_dofpt = Xh1->dof()->dofPointBegin();
    typename space_type::dof_type::dof_points_const_iterator en_dofpt = Xh1->dof()->dofPointEnd();

    for( ; it_dofpt != en_dofpt; ++it_dofpt   ) {
        kd_tree.addPointWithId( boost::get<0>( *it_dofpt ) , boost::get<1>( *it_dofpt ) , boost::get<2>( *it_dofpt ) ) ;
    }

*/
  //  typename MyKDTree::Element* M_tree;
    /*M_tree = */




/*
    std::cout<<"\n\n";
    std::cout<<kd_tree.points();
*/
}
/*
template<uint Dim,uint N>
void
Interpolate<Dim,N>::solve( sparse_matrix_ptrtype& D,
                  element_type& u,
                  vector_ptrtype& F )
{
    // solve the system, first create a vector U of the same size as
    // u, then call solve,
    vector_ptrtype U( M_backend->newVector( u.functionSpace() ) );
    // call solve, the second D is the matrix which will be used to
    // create the preconditionner
    M_backend->solve( D, D, U, F );
    // copy U in u
    u = *U;
} // Interpolate::solve


template<uint Dim,uint N>
void
Interpolate<Dim,N>::exportResults( element_type& U )
{
    Log() << "exportResults starts\n";

    exporter->step(0)->setMesh( U.functionSpace()->mesh() );

    exporter->step(0)->add( "u", U );

    exporter->save();
} // Interpolate::export
*/

// instantiation
template class Interpolate<2,2>;

